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ABSTRACT 

In the standard picture of cosmological structure formation, initially random-phase 



, fluctuations are amplified by non-linear gravitational instability to produce a final 

[t I ' distribution of mass which is highly non-Gaussian and has highly-coupled Fourier 

phases. We use the Zel'dovich approximation in one dimension to elucidate the onset of 
non-linearity including mode spawning, merging and coupling. We show that as grav- 
itational clustering proceeds, Fourier modes are spawned from parent ones, with their 
CO ' phases following a harmonic relationship with the wavenumbers. Spawned modes could 

^ , also merge leading to modulation of the amplitudes and phases which consequently 

■ breaks such harmonic relation. We also use simple toy models to demonstrate that bis- 

^■f^ , pectrum, Fourier transform of connected three-point correlation functions, measures 

' phase coupling at most at the second-order only when the special wavenumber-phase 

, harmonic relation holds. Phase information is therefore partly registered in bispectrum 

and it takes a complete hierarchy of polyspectra to fully characterize gravitational 
non-linearity. 

Q^. Key words: cosmology : theory - large-scale structure of the Universe - techniques: 



analytical 



1 INTRODUCTION 



Ij ' One of the main issues in cosmology is quantitative characterization of the large-scale structure of the Universe. The Universe 
, from present observations is highly non-linear on scales up to roughly 8/i~^Mpc. In the framework of the inflation paradigm, 
d ' the structure we see today is blown up from tiny inhomogeneities originated from quantum fluctuations. The cosmic microwave 
background has provided a strong evidence that the Universe was fairly homogeneous and isotropic with fluctuations one part 
in 10^ in the past. Such initial density fluctuations in the simplest inflationary universe constitute a Gaussian random field. 

Gaussian random Fields (Bardeen et al. 1986) are useful because of its analytical simplicity. One particularly interesting 
property of Gaussian random fields is that the real and imaginary part of the Fourier modes are both Gaussian distributed 
and mutually independent, or in other words, the Fourier phases are randomly distributed between and 27r. The statistical 
properties are then completely specified by its second-order statistics: its two-point correlation function ^(r) (with zero 
connected n-point correlation functions for n > 2), or alternatively, its power spectrum P{k). 

In the framework of gravitational instability, a perturbative method can be adopted at the early stage of density clustering 
(see e.g. Peebles 1980, Bernardeau et al. 2002). The linear perturbation theory is applicable when the density fluctuations are 
small i.e. its variance of density contrast {S^) < 1. The statistical distribution of the density fleld such as originally random 
phases remains invariant in the linear regime, except its variance increasing with time. 

The departure of the density fleld from the linear regime gives rise to phase coupling. Second-order statistics, such as 
power spectrum and two-point correlation function, throw away the fine details of the delicate pattern of cosmic structure. 
These details lie in the distribution of Fourier phases to which second-order statistics are blind. The evident shortcomings of 
P{k) can be partly ameliorated by defining higher-order quantities such as the bispectrum (Peebles 1980; Goroff, Grinstein, 
Rey & Wise 1986; Bernardeau 1992; Hivon et al. 1995; Matarrese, Verde & Heavens 1997; Scoccimarro 1997; Scoccimarro 
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et al. 1998; Scoccimarro, Couchman & Frieman 1999; Verde et al. 2000) or correlations of 5{k)'^ (Stirling & Peacock 1996). 
Higher-order correlations and polyspectra find this information term by term so that an infinite hierarchy is required for a 
complete statistical characterization of the fluctuation field. 

Due to the close connection between morphology and Fourier phases (Chiang 2001), there have been attempts in inves- 
tigating the behaviour of phases (Ryden & Gramann 1991; Soda & Suto 1992; Jain & Bertschinger 1998). They focus on 
phase shift from the original one. Scherrcr, Mclott and Shandarin (1991) develop a practical method on quantifying phase 
coupling. Chiang & Coles (2000) take phase difference between neighbouring modes on Shannon entropy to quantify phase 
information. Due to the circular nature of phases, a novel visualization method is also developed (Coles & Chiang 2000). 
Chiang, Coles & Naselsky (2002) and Chiang, Nasclsky & Coles (2002) develop a novel method, return mapping of phases to 
render phases onto a 'return map' to quantify phase associations between different Afc. Watts, Coles & Melott (2003) have 
shown that the probability distribution of phase difference between neighbouring modes displays a universal behaviour for 
clustering phenomenon. 

Since both Fourier phase coupling and non-zero higher-order correlation functions depict gravitational non-linearity. 
Watts & Coles (2002); Matsubara (Matsubara 2003) have shown the relationships between phase coupling and bispectrum 
which displays the lowest order non-linearity in density perturbations. Hikage, Matsubara & Suto (2003) have used phase sum 
from wavenumbers (fe vectors) forming a triangle as a descriptor of non-linear gravitational clustering. 

Despite the progress described above, efforts are yet to be made in understanding the direct relationship between gravi- 
tational clustering and phase coupling. In this paper we use onc-dimcnsional Zel'dovich approximation as a clustering scheme 
to investigate mode coupling when gravitational perturbations depart from the linear regime. This moment is related to the 
so-called quadratic density fields (Coles & Barrow 1987; Watts & Coles 2002 and references therein), in which the quadratic 
term induces non-linearity. The non-linearity caused by gravitational clustering can be described in two effects: mode spawning 
by the quadratic and higher-order terms and mode merging leading to modulation of the amplitudes and phases. These two 
effects are in action as cosmological density clustering proceeds. 

This paper is organized as follows. In Section 2 we re-visit some useful technical background: the commonly-used sta- 
tistical tools and theories: the power spectrum and two-point covariance functions, the linear theory of cosmological density 
perturbations, and the Zel'dovich approximation which in one dimension provides the insight on phase coupling. In Section 3 
we use such one-dimensional model as a clustering scheme for analyses of the onset of the following phenomena: mode spawn- 
ing and merging from gravitational clustering and their relationships with bispectrum. We further examine these phenomena 
by direct ID simulations into highly non-linear regime in Section 4. The discussions are in Section 5. For completeness in the 
Appendix we prove the exactness of the density evolution from the ID Zel'dovich approximation. 



2 TECHNICAL BACKGROUND 

2.1 Power spectrum and two-point covariance functions 

The mathematical description of an inhomogeneous Universe revolves around the dimensionless density contrast, S{x), which 
is obtained from the spatially-varying matter density p{x) via 

= (1) 
po 

where x is the comoving coordinate, po is the global mean density (Peebles 1980). The two-point covariance function, which 
measures the excess probability over Poisson distribution between a pair with separation r, is defined as 

ar) = {5{x)S{x + r)), (2) 

where the mean is taken over all points a;. It is also useful to expand the density contrast in Fourier series, in which 5 is 
treated as a superposition of plane waves: 

6{x) = S{k) exp(zfe • x). (3) 
The Fourier transform 5{k) is complex and therefore possesses both amplitude \S{k)\ and phase 0^ where 
5(fc) = |5(fe)|exp(»0fc). (4) 
The power spectrum is defined as 

(5(fei)5(fe2)) = (27r)'P(fc)5^(fei + fea), (5) 

which is the Fourier transform of the two-point covariance function via Wiener-Khintchin theorem. We can analogously define 
the bispectrum as the third-order moment in Fourier space: 

{5{ki)6{k2)6{k3)) = (27r)^B(fei, fe2, k3)5''{ki +k2+ fes), (6) 
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which is the Fourier transform of the connected three-point covariance function. 
2.2 The linear theory of density perturbations 

Gravitational instabilities is believed to be the driving force of large-scale structure formation. When the density perturbation 
is small, evolution of the density contrast can be obtaiucd analytically through the linear perturbation theory from 3 coupled 
partial differential equations. They are the linearized continuity equation, 

7n = — Vx-r, (7) 
at a 

the linearized Euler equation 
dv d 1 „ 1„ , 

-57 + -« = Va.p- -Va,0, (8) 

at a pa a 

and the linearized Poisson equation 
^(j) = 4:nGa^poS. 

In these equations, a is the expansion factor, p is the pressure, Vx denotes a derivative with respect to the comoving 
coordinates x, v — ax is the peculiar velocity and <j)(x,t) is the peculiar gravitational potential. Prom Eq.(7)-(9), and if one 
ignores pressure forces, it is easy to obtain an equation for the evolution of 5: 

6 + 2{'^)S -4nGpoS = 0. (10) 

For a spatially flat universe dominated by pressureless matter, po(t) = l/GirGt^ and Eq.(lO) admits two linearly independent 

power-law solutions 

S{x,t) = b±(t)SQ(x), (11) 

where do{x) is the initial density distribution, 6+(t) oc a(t) oc t^^^ is the growing mode and b-{t) oc is the decaying mode. 

As one can see from Eq.(ll), the growth depends only on time, so the growth of density distribution in the linear regime 
does not alter the phases. 



Va:^4> = 4nGa'^poS. (9) 



2.3 The Zel'dovich Approximation 

In order to examine analytically non-linear effects induced by gravitational clustering, we use the Zel'dovich approximation 
(1970) as a clustering scheme, which extrapolates the evolution of density perturbations into non- linear regime. This extrap- 
olation from the linear theory follows the perturbation in particle trajectories rather than in density fields. In the Zel'dovich 
approximation (ZA), a particle initially placed at the Lagrangian coordinate q is perturbed after a time t to an Eulerian 
coordinate x. The displacement of the particle simply depends on the constant velocity the particle has when it is kicked off 
its initial position, and can be written as b{t)u{q), so that 

r{q,t) = a{t)xiq,t) = a{t)[q + b{t)u{q)], (12) 

where r is the resultant physical coordinate, a{t) is the expansion factor and b{t) is the growing mode 6+(i) from the linear 
perturbation theory. According to this prescription, each particle moves with a constant velocity along a ballistic trajectory, 
which resembles Newtonian inertial motion. Note that the peculiar velocity according to the ZA is ax = a{t)b{t)u(q) . For an 
irrotational flow, u{q) can be expressed as a gradient of some velocity potential —V<to{q). 

We focus particularly on the applications of one-dimensional ZA. The ZA in ID provides an exact solution of density 
evolution (Buchert 1992) in that the evolution of planar collapse from the ZA has the same solution as from the Poisson 
equation until shell-crossing (Padmanabhan 1993). The ZA in ID is simplified to 

x{q,t) = q + b{t)u{q) = q-b{t)^^0^. (13) 
The density contrast can be derived from the conservation of mass pdx = podq : 



Po V ^1 



l-b{t)- 



dq^ 



dg2 

The velocity potential $0(5) can always be disintegrated as 



(14) 
(15) 
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^o{q) = '^AiCos{Xiq + ai). (16) 

i 

One can recover the solution to the linear theory from Eq.(14) by taking partial differentiation d/dt on both sides. Substituting 
(f^o/dq^ with b~^5/{d + 1), we reach 

+ A (17) 

Equation (17) provides us the insight into phase coupling from gravitational clustering. One important property of Eq.(17) 

is that the 5^ term is the "culprit" of the onset of non-linearity. Without the 5^ term we obtain the same solution as that 
from the linear theory, (5 oc 6 (oc t^^'^ if the Universe is matter-dominated), i.e. density fields grow only with time in the linear 
regime where the phases are unchanged (and uncorrelated if primordial Gaussianity is assumed). It is this quadratic term 6^ 
that the density field breaks away from the linear regime and produces gravitational non-linearity, hence phase coupling. In 
Eq.(17), the quadratic term is related to the "quadratic density field" (Coles & Barrow 1987; Watts & Coles 2002). 



3 MODE COUPLING AND BISPECTRUM 

In this section we use ID ZA to demonstrate spawning and merging of Fourier modes and phase coupling induced by 
gravitational clustering and its relation with bispectrum. For the analyses in this Section, we consider only an early stage of 
evolution, i.e. b{t) ^ 1 and assume that the curvature of the velocity potential is small such that x q and dx ~ dq. Fourier 
transform can then be performed in the Lagrangian coordinate q. Direct simulations of ID ZA for gravitational evolution is 
in Section 4. 



3.1 Mode spawning and mode merging 

Upon using ID ZA, we firstly choose as a toy model the velocity potential 



t 



^o{q) = -Ai cos(Aig -|- ai) - Ai cos(A2g -I- 02). (18) 
When 6(t) ^ 1, the density contrast can be expressed in terms of only the first order from Eq.(15) 

5 ~ (5*^^ = 6(t) ^^^^^^ = ai cos(Aig + ai) + 02Cos(A2g + 02) s 01 ^ ^ + 02 ^ ^, (19) 

where the round brackets hereafter denote cosine functions, a\ = 6(i)AiAf and 02 = h(t)A2y^. After Fourier transform in 
Lagrangian coordinate we have only two modes with wavenumbers Ai and A2, and phases ai and a2, respectively. So these 
are the only Fourier modes the moment the clustering process begins (which we call the "parent modes" ) . It is clear that the 
first order does not display phase coupling. When the second-order term becomes comparable. 



5 . 5«+5(^) = 6(t)(^)+6^(t) 



dg2 



ci-i+0,2, f \ , f ^2 \ , 2 ( 2Ai \ 2 ( 2A2 A / Ai + A2 \ / Ai — A2 



+ ai ( 








+ ^2 f 






, "1 / 




, ^2 / 



+ ai[ \+ai[ +aia2 ^' " ) + a.a^ [ '"^ . (20) 

\ 2ai J \ 2a2 J \ + 02 J \ ai - 012 J 

Since the density contrgist in ID ZA can be expressed as a power series of d^^o/dq^ in Eq. (15), the second-order term 
reflects the quadratic density fields (Coles & Barrow 1987; Watts & Coles 2002). This second-order term spawns modes with 
wavenumbers that are from combinations of any 2 wavenumbers from the parent ones, i.e. Fourier modes with wavenumbers 
2Ai, 2A2, Ai -|- A2 and Ai — A2 spawned from combination of Ai and A2. One important feature of quadratic density fields is 
that the phases of the spawned modes follow the same kind of harmonic relationship as the spawned wavenumbers, which we 
call hereafter "wavenumbor-phaso harmonic relation". Such relation subsequently forms phase associations between Fourier 
modes and is crucial for bispectrum analysis. 

We can generalize the velocity potential in ID ZA as a sum of cosine functions $0 = — Ai co8{Xiq + ai). The density 
contrast of the first order and that up to the second order now become, respectively, 

S^^^ = ^ ai cos(Ai +ai) = Y,'^4 ) . (21) 
where o, = b{t)AiXi, and 

^ We put minus sign for both cosine functions in the velocity potential to make the analytical form neat, assuming both Ai and A2 
positive, otherwise the phases would have a shift by tt. 
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i ^ ^ jk "-^ ' 



(22) 



If primordial Gaussianity is assumed, the phases a, of the initial velocity potential are uniformly random and independently 
distributed between and 27r. Below we categorize the effects induced by quadratic density fields. 

• Mode spawning : In Eq.(22) beyond the first order, new modes are spawned from the 5'^^ and correlated phases 
are created following the spawned modes. The wavenumbers of the spawned modes are formed from combination of any 2 
wavenumbers of the first order (the same as the toy model) and the phases follow the wavenumber-phase harmonic relation. 
Some terms may appear oven at the earliest stage (when b{t)AjAk{Xj\k/\i)'^ > 2Ai). Those are terms usually involving 
high-frequency modes, i.e. high or A^. Such phase coupling can manifest itself through phase mapping (Chiang, Coles & 
Naselsky 2002). For example, a short sequence of modes is formed with a constant difference in wavenumber Ak = \k' 



Ojak 



A. 



ak 



,aj 



ajttk ( Aj + Afc 



(23) 



where the middle mode is taken from the cj'^-*. Due to the wavenumber-phase harmonic relation, this phase sequence has a 
constant phase difference Acj) = Ok and can be mapped along a line parallel to the diagonal through phase mapping technique. 
Such coupling between modes with large Afe is discussed in Chiang, Coles & Naselsky (2002). 

• Mode merging : The newly-spawned modes from the second order are not all independent but could merge with modes 
of the same wavenumbers in the first order. Take the following two modes from Eq.(22) as an example, one from the j'^' and 
the other from the 5^^^: 



Xi 



(ijak 
2 



A,- 



Afc 

ak 



(24) 



If A, 



Afc = Ai, i.e. the wavenumber of the newly-spawned mode coincides with the one in the first order, and when b{t) 

ajak/2, these two modes merge to form a new mode: 



reaches a stage b{t) 

2ai cos[(aj — Ok — ai)/2] 



2{Ai/AjAk){Xi/XjXk)^ such that ai 

A, 



(25) 



where the amplitude is modulated and the phase is shifted. Such modulation in amplitudes and phases proceeds gradually 
and the consequence is that the wavenumber-phase harmonic relation is broken. Thus not all spawned modes will enjoy 
the harmonic relation. Note that in this case the evolution can still be at very early stage b{t) <§; 1 as long as XjXk/Xi » 
{Ai/AjAky^'. 



We can also look at the density contrast up to the third order 

A, 
a 



(2) 



?(3) 



+ 



+ 



^ai 

i 

Eatam,an 
4 



2 

jk 

Xi -|- Am + A, 
Oit + am + a, 



[( 


' A, 


— Afe 


H 


Xj + Xk \ 






— «fe 




, aj +ak J 



Xi + A„ 



Ar, 



ae + Um — a„ 



Xt — Xn 

ae — an 



A„ 
+ an 



Xt 
at • 



■ Aj7 

■ an 



Xji 
a„ 



(26) 



Once again new modes are spawned, more complicated phase coupling appears and the spawned modes also merge with parent 
modes. We shall look at mode spawning and coupling at higher-order modes through direct simulations. 



3.2 Bispectrum and quadratic phase coupling 

Bispectrum is the lowest-order statistic sensitive to the structure generated by gravitational clustering (Peebles 1980). It is 
defined as two-dimensional Fourier transform of the connected three-point auto-covariance function C, 

C(ri,'r2) = {5{x)5(x + ri)5{x + T2)). (27) 

The bispectrum and higher-order polyspectra vanish for Gaussian random fields, but in a non-Gaussian field they may be 
non-zero. The usefulness of these and related quantities therefore lies in the fact that they encode some information about 
non-linearity and non-Gaussianity. 

A certain form of phase relationship produces non-zero bispectrum. To see this, we simpfify the of Eq.(20) by taking 
only the 2nd, 3rd and 6th term and denote it as the case I (Coles & Chiang 2001), 

5i = aif \+aJ )+a^aJ + V (28) 



a\ I \ a2 I V Qi -|- 012 
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For comparison we introduce the case II: 

On = 01 + + oia2 , (29) 



ai J \ 012 J \ 0:3 

where ai, a2 and as are random. The case I displays quadratic phase coupling whereas the case II exhibits no phase association. 
One can see that {5i) = {5u) = 0, and it is straightforward to show that the covariances for both cases are equal, 

2 2 22 

6(r) = {5i{x)5i{x + r)) = 6i(r) = y cos(Air) + ^ cos(A2r) + ^ cos[(Ai + Aajr], (30) 

so are their power spectra. The above demonstrates that second-order statistics are blind to any phase associations. 
The reduced three-point covariance function for the case I is 
^2^2 

Ci(n,r2) = ^-^[cos(A2ri + Air2) + cos(Airi -I- A2ri - Air2) -|- cos(Airi + A2ri - A2r2) 

-|-cos(Airi + Mvi) + cos(Air2 + \2r2 — Airi) -|- cos(Air2 + \2r2 — A2ri)], (31) 
and for arbitrary ai , a2 and as we get 

Cii(r-i,r2) =0. (32) 

One can see that the spawned modes of the second order which benefit from the wavenumbor-phasc harmonic relation induce 
non-zero The bispectrum, B{ki,k2), the two-dimensional Fourier transform of for case II is Bu{ki,k2) = trivially, 
whereas -Bi(fci, fc2) consists of a single spike located somewhere in the region of {ki, ^2) space defined by k2 > 0, fci > k2 and 
ki + k2 < TT. If Ai > A2 then the spike appears at ki = Ai, k2 = A2. Thus the bispectrum is the lowest order of polyspectra 
that measures the phase coupling induced by quadratic non-lineaxities when the wavenumber-phase harmonic relation holds. 
However, bispectrum is trivially zero again if such harmonic relation is broken. Note that this relation is easily broken by 
mode merging when a third phase is involved. 

Following the same lino of thought, we can extend this toy model from Eq.(26) to the case III, 

5iu = ai[ +0,2 + "3 "I A — _L _L ' (33) 

\ ai J \ a2 J \ a-i J 4 \ ai + a2 + as J 

and the case IV for comparison, 

5,v = a,(^^J+a2^^J+a3(^^J+^(^ )' ^'"^ 

where ai, 02, Q3 and 04 are random. The bispectrum for both cases are zero, but the trispectrum, the three-dimensional 
Fourier transform of the reduced four-point covariance function, can pick up cubic phase coupling for case III (24 terms), 
whereas it is again trivially zero for case IV. The third-order term produces cubic phase coupling to which bispectrum is blind. 

Therefore, bispectrum only measures phase coupling at the second-order. During density clustering, however, mode 
spawning will produce the bispectrum from case I, Ci(''ii''2) is not always non-zero. Mode merging breaking the wavenumber- 
phase harmonic relation results in zero bispectrum. The same would happen to all polyspectra. 

Although we only demonstrate with ID toy model the relationship between bispectrum and quadratic phase coupling, a 
full treatment of such relationship is elaborated in Watts & Coles (2002): a non-zero bispectrum can only be produced when 
the k vectors form a triangle in fe-space, wherein the same (quadratic) wavenumber-phase harmonic relation has to hold. 
Phase coupling is not registered only in bispectrum, but in all polyspectra. 



4 SIMULATIONS 

In this Section, we perform ID ZA numerically in order to demonstrate mode spawning and coupling into highly non-linear 
regime. Due to limitation on integration in Lagrangian coordinate, in the previous Section we can analyze only the onset of 
non-linearity by assuming b{t) <C 1. In the following simulations we assume the periodic boundary condition, 

x{q + 2jTv) =x{q). (35) 
4.1 Non-interactive and interactive pcirent modes in the simulations 

Here we supply three different initial velocity potential functions to the simulations. We choose the velocity potential field 
with the following form for Fig. 1 and 2: 

$0 = -cos(7g-|-6») -cos(ll9-|-</>), (36) 
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Figure 1. The amplitude and phase evolution of ID ZA simulation with $o = — cos(7q) — cos(llq + 7r/16) for wavenumber fc = 1 — 49. 
We specifically designate 7 small panels in one row so that in each column the increment of the wavenumber is 7. Ea<;h time step is 
represented by a diamond symbol. In each panel (for one Fourier mode) the amplitude (xlO^) is shown in logarithmic scale by the 

distance to the origin, and the phase by the angle against the positive x axis. Thus the dotted circle with radius 5 units indicates the 
level of amplitudes with 10^^. The radial dotted lines arc angles with increment of 7r/16 to show that phases align with angles of multiple 
of 7r/16. The dash lines arc the predicted phases according to Table 1: the wavcimmber-phase harmonic relation. One can see the phases 
of spawned modes follow this harmonic relation because they have only one primary combination from parent modes, i.e. non-interactive. 



where B and are and 7r/16 in Fig. 1 and 7r/4 and 7r/7 in Fig. 2, respectively. As described in previous Sections, once 
the clustering proceeds, Q and will feed as phases of two parent modes with wavenumber 7 and 11 in these cases: a^ and 
Qii. These specially chosen wavenumbers Ai = 7 and A2 = 11 for both Fig. 1 and 2 are "non-interactive" in the sense 
that the phases of spawned modes have only one primary combination from the parent ones. If the parent modes are "non- 
interactive" , the phases will eventually follow the harmonic relationship (see Table 1). So the phase of the spawned mode 
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with, for example, wavenumber fc = 25 ( = 2x7 + 11) can only come from addition of the phases of the parent modes, in 
other words, a2s = 2 x ar + an. 

In addition to the two simulations of "non-interactive" mode spawning, in Fig. 3 we also show the simulation with one 

extra "interactive" mode: 

$0 = - cos(7g) - cos(llg + n/16) - cos(14g + n/4). (37) 

Figure 3 with the extra mode A3 = 14 = 2Ai illustrates parent-mode interactions leading to change of phases during the 
clustering process. One can see now, with this extra parent mode, the spawned wavenumber k = 25, for example, can have 
two combinations: 2 x 7 + 11 and 11 + 14. So the phase of this spawned mode 025 can now come from either 2a7 + an or 
an + ai4. 

We using the following representation in order to display simultaneously the amplitude and phase evolution for each 
Fourier mode. Each small panel represents one Fourier mode with denoted wavenumber. Wc divide the simulations (up to the 
shell crossing) into 20 time steps with each step represented by one diamond symbol. In each panel the amplitudes of different 
time steps are shown in logarithmic scale by the distances to the origin, and the phases by the angles against the positive x 
axis (similar to a complex-plane representation but only the magnitudes in logarithmic scale). We have to multiply all the 
amplitudes by 10® so that they are always positive after taking the logarithm^ . Thus the dotted circles with radius 5 units in 
each panel indicate a level of Fourier amplitudes 10~*. 

In Fig. 1 the radial dotted lines arc angles with increments of n/16 to show that phases aligned with angles of multiple 
n/16 (which is due to combination of 07 = and an = n/16). The dash lines are the predicted phases from wavenumber-phase 
harmonic relation according to Table 1. 

One can easily see in Fig.l that the parent modes k = 7 and 11 have higher initial amplitudes than other modes for 
the start. We specifically designate 7 small panels in one row so that in each column the increment of the wavenumber is 
7. Because wavenumber 7 is one of the parent modes, phases of the spawned modes are formed by addition of ay, following 
wavenumber-phase harmonic relation. For instance, for the wavonumbers of the Fourier modes 18 {— 11 + 7), 25 (= 11 + 2 x 7), 
32 (= 11 + 3 X 7), 39 (= 11 + 4 X 7), 46 (= 11 + 5 x 7), their phases are formed from an + ar, an + 2 x 07, an + 3 x a7, 
an + 4 X a7 and an + 5 x a7, an increment of a7 (see Table 1 for the details of phase combination from the parent modes). It 
is therefore easy to see that the multiples of wavenumber 7: 7, 14, 28, 35, 42 have the same phase, as well as 8, 15, 22, 29, 36, 
43. This is due to the fact that the phases formed by the increment of a7 is unchanged when a7 = in Fig.l. Phases of most 
modes eventually align with multiples of 7r/16 because of the following reasons: they follow the wavenumber-phase harmonic 
relation and the phase of one parent mode is 0. Modes with wavenumbers A; = 1, 2, 5, 6, 9 and 13 do not align with multiples 
of 7r/16 and their amplitudes cannot reach 10"''. 



4.2 Analytical account for the simulations 

To account for the interaction of the parent modes leading to mode merging, we can start from Eq. (14). We have 

J —TV \ J —TK J —TV 



(38) 



where u{q) = —d<^o{q) / dq. We can give the same treatment to Eq. (38) as Eq. (15) by expressing it in terms of a power series 
of 6(t): 



p=0 



(-ikbY 1 
pi 2^ 



dq e-''"'u{qy. 



Following Soda & Suto (1992), we can denote the Fourier transform of u{q) as Uk 
dq u{q)e~'^'"^, 



1 r 

= 2^1 ' 

J —TV 



so the Fourier modes 6k can be expressed as 

{-tkbf 



5k = -ikbuk + Uk^Uk^S^lk - (fci + k2)] + ^ ^ UkiUk2Uk3&° [k - {ki_ + k2 + k^)] + 



(39) 



(40) 



(41) 



where (5° denotes Dirac-delta function. For u{q) = —d^o/dq = —'^jAj\j8in{\jq + aj) from the velocity potential $0 = 
-Sj^jCos(Ajg + aj), 

■t This is to avoid the degeneracy between negative amplitudes (after taking the logarithm) with positive phases and positive ampUtudes 
with opposite phases. 
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-3ar + 2aii + w 



5aj — 3aii + tt 



03 = 
2a7 — ail 



a4 = 
-ar + ail 



as 



ae • 



-Aa-j + 3aii 



-4a7 + 2aii 



as = 
-2a7 + 2aii 



ag 



6a7 — Sail + TT 



aio = 
3a7 — ail 



ail 



ai2 = 
— 3a7 + 3aii 



ai3 : 



5a7 — 2aii 



ai4 = 
2a7 



-a7 + 2aii 



ai6 = 
-4a7 + 4oii 



ai7 = 
4a7 — ail + TT 



ai8 = 
07 + ail 



ai9 = 
-2a7 + 3aii 



a-za = CX21 = 

6a7 — 2aii 3a7 



"22 = 
2<ui 



"23 = 

-:5m7 + Inii 



"24 = 
■^)"7 - nil + TT 



"25 — 

2"7 + nil 



"26 = 
"7 + •■-!<Ui - : 



"27 = "28 = 

— in 7 + on 11 in 7 



"29 = 
"7 + 2aii 



"30 = 

-2a7 + 4aii 



"31 = 
6"7 — "11 + TT 



"32 = 

3a7 + ail 



0:33 = 
3aii 



"34 = a35 = 

-3o7 + 5aii + TT 5"7 



"36 = 

2a7 + 2aii 



"37 = 

-ar + 4aii + tt 



"38 



— 4a7 + 6aii 



a3g = a4o = a4i = 042 = 

4a7 + ail a7 + 3aii — 2a7 + 5aii 6a7 



"43 = 
3a7 + 2aii 



a44 = 
4aii 



"45 = 
-3a7 + 6aii + n 



a46 = 
5a7 + ail 



a47 = 
2a7 + 3aii 



a48 = "49 = 

-ar + 5aii + tt | 7a7 



Table 1. The predicted phases from the hariiiouic relation for parent wavemiinljer 7 and 11 for Fig.l and 2. The terms with odd 
number of complex conjugate give rise to an extra phase shift by tt, apart from the harmonic relation. The spawned modes that 
have combinations from the parent wavenumbers (with lower orders) will grow, with phases following the same harmonic relation. Those 
suppressed are modes which need high orders of combinations from wavenumbers of the parent modes (those in frames, e.g. fe = 13 and 
38). 



uk = V A, A,e»°^ <5^(A, - k) = -ii^e'«'= = (42) 
2i ^ — ' 2i 2 

j 

and M_fe is simply its complex conjugate. 

Equation (41) and (42) provide us the insight into mode spawning and mode merging in gravitational clustering. The 
Dirac-delta functions dictate the spawning of Fourier modes, both the amplitudes and phases. If, for any Fourier mode k, there 
is only one combination from the parent modes, its phase shall not change and the amplitude grow steadily with b(t)" 
(according to the order of this combination). If, on the other hand, there is more than one set of combination for wavenumber 
k [cf Eq. (24)], they will interact so that the phase ak will be changed accordingly. 

To be more specific, in our case of two cosine functions in <I>()(g) (in Fig.l and 2) from which only and ux^ exist 
(Af, > Ao), we can write down explicitly what constitutes the Fourier mode of wavenumber k: the first term in Eq. (41) is 
simply AkXkkbe^"'' /2. The second term has all the combinations of any two from Aa and Aj in the Dirac-delta functions: 2Aa, 
2Ai,, Aa + Ai, and Aa — Xb so that 



{-ikbf 
2! 

~ 2 

+ ^L^e'(^«^a+^><5^(fe - 2Aa) + ;l^e'(2"^6+-)j«(fc _ 2A,) 



UkiUk^S'^lk - {ki + fe)] 

fci ,fc2 

^Aq^AfcAgA!) _z(a;,„+a;,^+7r) jc-D/i, a a.\ , ^A„^A(,AaA6 



-e'("^»+"^i>+"^5^(fc - Aa - Xb) + ' ;'-^°^'' e'<"^6-"A,)jg(fc + Xa- Xb) 



(43) 



If the $0 has more than two cosine functions, the second term of Eq. (41) then shall have all the combinations of any 2 from 
$0- The third term involves combinations of any 3 from Aa and At, i.e. 3Aa, 3Ai,, 2Aa ± At, Aa ± 2Ai,. It is also worth noting 
that unless the higher-order terms exist (i.e. — {ki + A;2 + ^3)] = 1), the phases of the Fourier modes are fixed and the 
amplitudes grow with 

For the interactive parent modes in Eq.(37) (Fig. 3), where A3 = 2Ai, the second term of Eq.(41) has one Dirac-delta 
function involving A3-I-A2 whereas in the third term there is another involving 2A1-I-A2. The resultant phase of this wavenumber 
depends on the amplitudes of the terms involving the two Dirac-delta functions. 

Note that for wavenumber A: = Ai, — Aa of Eq.(43) the phase is qa^ — -\- n. This extra phase shift by tt is due to odd 
number of complex conjugate of Uk in Eq. (41). 

In Table 1 wo show the predicted phases of the spawned Fourier modes from the parent modes for Fig. 1 and 2. The 
wavenumbers of the parent modes are 7 and 11. The spawned modes that have the combination of the parent wavenumber 
(with lower orders) will grow following the wavenumber-phase harmonic relation. Those which are suppressed (or have not 
grown before shell crossing) are modes which need high orders of combinations in wavenumbers or which have no exact 
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Figure 2. The amplitude and phase evolution (wavenumber k = 1 — 49) of ID ZA simulation with ^>o = — cos(7q + 7r/4) — cos(llg+7r/7). 
In each panel the amplitude (xlO^) is shown in logarithmic scale by the distance of ea<;h symbol to the origin, and the phase by the angle 
against the positive x axis. Thus the dotted circle in each panel indicates amplitude with 10~*. The dash lines indicate the predicted 
phases according to Table 1: the wavenumber-phase harmonic relation. 



combination. The phase of wavenumber 24: a24, for example, can have combinations of Sa? — Qii or —Gaj + Qaii . The former 
is of the 6th order and the latter is the 12th order in terms of h{t) series in Eq.(41). The low-order combination of wavenumbers 
prevails at the early stage and if it is dominant in amplitude, the other combinations will have little influence. When the 
parent modes are "interactive", the phase sets out with a value from lower-order combination. In Fig. 3 we have 3 parent 
modes Ai = 7, A2 = 11 and A3 = 14. The phase of mode fe = 18 is spawned from ar-l-aii in the beginning and Sar + an —an 
at later stage. It is only at very later stage of evolution that the 5th order can catch up with the 2nd order (see Fig. 3). 

Having such rules in mind, one can clearly see why the parent low-frequency modes (and their phases) have decisive 
influence on the mode coupling. As shown in Chiang & Coles (2000), convergence of phase difference Dk = ipk+i — 't>k at high 
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Figure 3. The amplitude and pliase evolution of ID ZA simulation with $o = — cos(7q) — cos(llg + 7r/16) — cos(14g + 7r/4). As in Fig. 1 
the radial dotted lines are angles with an increment of 7r/16 to show that phases align with angles of multiple of 7r/16. One can see now, 
with the extra parent mode A3 = 14 interacting, the phases are shifted considerably. Note also that in this case, the wavenumber-phase 
harmonic relation is broken. 



fc is related to the location of the highest density peak. The wavenumbers of parent lowest- frequency mode (fc = 1) can easily 
be added (or subtracted) to form precisely the spawned wavenumbers, particularly for high-frequency modes, so the phase of 
the lowest k mode would inevitably take part in most of the hight fe mode spawning. Hence the increment in phases (i.e. D}2) 
along k axis would be the phase from the lowest frequency mode (fe = 1). 

Note also that for more generic situations where all modes in $0 have non-zero amplitudes, they are definitely "interactive" , 
in particular with the presence of fundamental mode fc = 1. Therefore, the evolution of amplitudes and phases shall have a 
complicated picture with mode merging with changing of Fourier amplitudes and phases. 
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5 DISCUSSIONS 

In this paper, we use ID ZA to illustrate some non-linear effects induced by gravitational clustering: mode spawning and mode 
coupling. We first demonstrate the onset of non-linearity by expanding the density contrast in terms of power series of b{t). The 
second-order of such expansion reflects quadratic density fields, from which modes are spawned with phases following a special 
wavenumber-phase harmonic relation. We have also illustrated with toy models that the widely-used statistic, bispectrum, can 
only pick up the quadratic phase coupling and is blind to cubic and higher-order phase coupling. We further perform direct 
ID ZA simulations in order to demonstrate those effects entering non-linear regime. The wavenumber-phase harmonic relation 
holds as long as the wavenumbers of parent modes are not "interactive" , as shown in Fig. 1 and 2. The complexity of Fourier 
amplitude and phase evolution comes from the continuous interactions between modes and their merging and coupling with 
the spawned modes. 

Recently it is also claimed that bispectrum measures phase correlations (Komatsu et al. 2003), which is partially correct. 
We have demonstrated with toy models that bispectrum measures only quadratic phase coupling under the condition that the 
wavenumber-phase harmonic relation holds, as pointed out in details in Watts & Coles (2002). Upon using bispectrum as a 
test of non-Gaussianity, a quadratic field such as the form 5{x) + ed^{x) could give a zero bispectrum due to the interaction 
between the spawned modes in ed'^{x) with the parent ones in S{x), hence resulting in a false signature of Gaussianity. 
Moreover, the blindness of bispectrum to higher-order phase coupling indicates that a complete hierarchy of polyspectra is 
required to fully characterize the statistical properties of a fiuctuation field. 

The initial power spectral index decides the clustering morphology, which is closely related to Fourier phases (Ghiang 
2001). It is therefore interesting to examine the link between Gaussian random fields with different initial power spectral 
index, the only available information, and evolved morphology with information encoded in phases rather than in amplitudes. 
The ID ZA simulations in this paper can be considered as from initial white noise power spectrum with spectral index n = 
(as the coefficients of cosine functions in the velocity potential arc all unity) . It is clear that mode spawning is independent of 
amplitudes of the parent modes, but mode merging and coupling is rather dependent on the amplitudes of spawned modes, 
which are originated from parent ones. How the information is transferred from Fourier amplitudes to different phase coupling 
configurations will be examined in the next paper. 
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